Day10:オープンデータの可視化
課題の回答例
【課題】限界集落の可視化
e-Stat(統計GIS)の平成22年国勢調査(少地域)から適当な地域を選び（出身地でも国立市でも、鯖江市以外なら良い）、「地域別高齢者(65歳以上)人口比率」がみえるように可視化してください。
可視化例：
if(!require(maptools)){
install.packages("maptools")
require(maptools)
}
if(!require(ggmap)){
install.packages("ggmap")
require(ggmap)
}
# ggplot2はggmap 読み込み時に自動的に読み込まれるが念のため
if(!require(ggplot2)){
install.packages("ggplot2")
require(ggplot2)
}
#日本語フォント設定（Win/Mac両対応)
if(.Platform$OS.type=="windows")
windowsFonts(yugo=windowsFont("Yu Gothic"))
if(capabilities("aqua"))
quartzFonts(yugo=quartzFont(rep("YuGo-Medium",4)))
#福井県鯖江市、平成22年国勢調査(少地域)2010/10/1
# shapefileの読み込み
shp<-maptools::readShapePoly("A002005212010DDSWC18207/h22ka18207.shp")
# data.frameへの変換
dat <- broom::tidy(shp)
# 行番号を新規列項目として追加
shp@data$id <- rownames(shp@data)
# 「dat」と「shp@data」とを列「id」をキーに結合
shp.df <- dplyr::inner_join(dat, shp@data,by = "id")
# 列「KEY_CODE」を数値に変換
shp.df$KEY_CODE<-as.numeric(as.character(shp.df$KEY_CODE))
# 統計データの読み込み
csv<-read.csv("tblT000573C18207.txt",stringsAsFactors = FALSE,fileEncoding = "SJIS")
# 内容確認
head(csv)## KEY_CODE HYOSYO CITYNAME NAME HTKSYORI HTKSAKI GASSAN
## 1 NA NA NA NA
## 2 1.8207e+04 1 鯖江市 NA NA
## 3 1.8207e+08 2 鯖江市 本町 NA NA
## 4 1.8207e+10 3 鯖江市 本町１丁目 NA NA
## 5 1.8207e+10 3 鯖江市 本町２丁目 NA NA
## 6 1.8207e+10 3 鯖江市 本町３丁目 NA NA
## T000573001 T000573002 T000573003 T000573004
## 1 総数、年齢「不詳」含む 総数０〜４歳 総数５〜９歳 総数１０〜１４歳
## 2 67450 3344 3604 3595
## 3 871 26 35 36
## 4 184 10 6 10
## 5 218 2 7 13
## 6 121 5 4 6
## T000573005 T000573006 T000573007 T000573008
## 1 総数１５〜１９歳 総数２０〜２４歳 総数２５〜２９歳 総数３０〜３４歳
## 2 3335 2768 3603 4487
## 3 30 20 26 50
## 4 4 5 10 11
## 5 10 7 7 7
## 6 2 4 3 16
## T000573009 T000573010 T000573011 T000573012
## 1 総数３５〜３９歳 総数４０〜４４歳 総数４５〜４９歳 総数５０〜５４歳
## 2 5303 4291 4038 3811
## 3 45 30 36 57
## 4 10 7 4 10
## 5 13 7 11 12
## 6 6 3 7 11
## T000573013 T000573014 T000573015 T000573016
## 1 総数５５〜５９歳 総数６０〜６４歳 総数６５〜６９歳 総数７０〜７４歳
## 2 4421 5373 4126 3470
## 3 52 81 65 63
## 4 14 20 16 13
## 5 11 12 17 17
## 6 8 14 8 7
## T000573017 T000573018 T000573019 T000573020
## 1 総数１５歳未満 総数１５〜６４歳 総数６５歳以上 総数７５歳以上
## 2 10543 41430 15371 7775
## 3 97 427 347 219
## 4 26 95 63 34
## 5 22 97 99 65
## 6 15 74 32 17
## T000573021 T000573022 T000573023 T000573024
## 1 男の総数、年齢「不詳」含む 男０〜４歳 男５〜９歳 男１０〜１４歳
## 2 32507 1686 1830 1846
## 3 386 14 21 14
## 4 84 4 2 4
## 5 93 1 5 5
## 6 53 3 1 1
## T000573025 T000573026 T000573027 T000573028
## 1 男１５〜１９歳 男２０〜２４歳 男２５〜２９歳 男３０〜３４歳
## 2 1768 1292 1761 2211
## 3 15 9 14 25
## 4 3 2 6 6
## 5 5 2 4 3
## 6 - 3 1 9
## T000573029 T000573030 T000573031 T000573032
## 1 男３５〜３９歳 男４０〜４４歳 男４５〜４９歳 男５０〜５４歳
## 2 2723 2141 1992 1867
## 3 22 14 17 24
## 4 7 1 4 1
## 5 6 4 4 6
## 6 2 2 3 6
## T000573033 T000573034 T000573035 T000573036 T000573037
## 1 男５５〜５９歳 男６０〜６４歳 男６５〜６９歳 男７０〜７４歳 男１５歳未満
## 2 2124 2595 2090 1598 5362
## 3 24 37 35 27 49
## 4 6 11 9 6 10
## 5 6 4 10 7 11
## 6 2 6 3 4 5
## T000573038 T000573039 T000573040 T000573041
## 1 男１５〜６４歳 男６５歳以上 男７５歳以上 女の総数、年齢「不詳」含む
## 2 20474 6597 2909 34943
## 3 201 136 74 485
## 4 47 27 12 100
## 5 44 38 21 125
## 6 34 14 7 68
## T000573042 T000573043 T000573044 T000573045 T000573046
## 1 女０〜４歳 女５〜９歳 女１０〜１４歳 女１５〜１９歳 女２０〜２４歳
## 2 1658 1774 1749 1567 1476
## 3 12 14 22 15 11
## 4 6 4 6 1 3
## 5 1 2 8 5 5
## 6 2 3 5 2 1
## T000573047 T000573048 T000573049 T000573050
## 1 女２５〜２９歳 女３０〜３４歳 女３５〜３９歳 女４０〜４４歳
## 2 1842 2276 2580 2150
## 3 12 25 23 16
## 4 4 5 3 6
## 5 3 4 7 3
## 6 2 7 4 1
## T000573051 T000573052 T000573053 T000573054
## 1 女４５〜４９歳 女５０〜５４歳 女５５〜５９歳 女６０〜６４歳
## 2 2046 1944 2297 2778
## 3 19 33 28 44
## 4 - 9 8 9
## 5 7 6 5 8
## 6 4 5 6 8
## T000573055 T000573056 T000573057 T000573058 T000573059
## 1 女６５〜６９歳 女７０〜７４歳 女１５歳未満 女１５〜６４歳 女６５歳以上
## 2 2036 1872 5181 20956 8774
## 3 30 36 48 226 211
## 4 7 7 16 48 36
## 5 7 10 11 53 61
## 6 5 3 10 40 18
## T000573060
## 1 女７５歳以上
## 2 4866
## 3 145
## 4 22
## 5 44
## 6 10# 先頭行の削除
csv <- csv[-1,]
dat <- data.frame(KEY_CODE=csv$KEY_CODE,TOTAL=csv$T000573019)
shp.df2 <- dplyr::inner_join(shp.df,dat,by="KEY_CODE")
shp.df2$TOTAL <- as.numeric(as.character(shp.df2$TOTAL))
shp.df2$RATE <- shp.df2$TOTAL/shp.df2$JINKO
# Inf(無限)となっていたものを0に置換(0除算結果の修正)
shp.df2$RATE[is.infinite(shp.df2$RATE)] <- 0
#NA(欠損値)を0に置換
shp.df2$RATE[is.na(shp.df2$RATE)] <- 0
#ラベル表示の準備
#文字コード変換(SJIS > UTF8).Windows環境の場合は不要．
shp.df2$MOJI <- iconv(as.character(shp.df2$MOJI),from="SJIS",to="UTF8")
#ラベル座標とラベルの設定
shp.text <- data.frame(shp.df2$X_CODE,shp.df2$Y_CODE,shp.df2$MOJI)
#重複を除去する
shp.text <- unique(shp.text)
names(shp.text) <- c("lon","lat","name")
shp.text$name <- as.character(shp.text$name)
#GoogleMapを取得する
map <- ggmap::get_map(location=c(136.2302503,35.9633669),zoom=12,maptype="roadmap")
#GoogleMapを描画する（ggplot2の初期化も実施している）
gm <- ggmap::ggmap(map)
#市町村の領域を描画.colorで境界線の色を指定．alphaで透明度を指定
g <- gm + ggplot2::geom_polygon(data = shp.df2, aes(x = long, y = lat, group = id,fill=RATE),color="gray80",alpha=0.7)
#塗り色を変更する
g <- g+ ggplot2::scale_fill_gradientn(colours=rev(heat.colors(100)))
#タイトルを指定する．
g <- g +ggplot2::ggtitle("福井県鯖江市の高齢者比率")+ggplot2::theme_gray(base_family = "yugo")
plot(g)
【発展課題】市営駐車場の可視化
以下の3つのデータを組合せた可視化にトライして下さい．
- 鯖江市のshapeファイル，
- 「データシティさばえ」の市営駐車場データ(http://www3.city.sabae.fukui.jp/xml/parking/parking.xml)
- 点の大きさで駐車場の最大収容数を表現してください
- Google Map
完成イメージは以下のとおりです。
if(!require(maptools)){
install.packages("maptools")
require(maptools)
}
if(!require(ggmap)){
install.packages("ggmap")
require(ggmap)
}
# ggplot2はggmap 読み込み時に自動的に読み込まれるが念のため
if(!require(ggplot2)){
install.packages("ggplot2")
require(ggplot2)
}
#XML処理に必要
if(!require(XML)){
install.packages("XML")
require(XML)
}
#shapefileを読み込む
shp<-maptools::readShapePoly("A002005212010DDSWC18207/h22ka18207.shp")
#data.frameに変換
dat <- broom::tidy(shp)
# 行番号を新規列項目として追加
shp@data$id <- rownames(shp@data)
# 「dat」と「shp@data」とを列「id」をキーに結合
shp.df <- dplyr::inner_join(dat, shp@data,by = "id")
# 列「KEY_CODE」を数値に変換
shp.df$KEY_CODE<-as.numeric(as.character(shp.df$KEY_CODE))
# 統計データの読み込み
csv<-read.csv("tblT000573C18207.txt",stringsAsFactors = FALSE,fileEncoding = "SJIS")
csv <- csv[-1,]
# 「csv$KEY_CODE,csv$T000573019」からなるdata.frameを作る
dat <- data.frame(KEY_CODE=csv$KEY_CODE,TOTAL=csv$T000573019)
# 「dat」と「shp.df」とを列「KEY_CODE」をキーに結合
shp.df2 <- dplyr::inner_join(shp.df,dat,by="KEY_CODE")
# 数値に変換（計算のため）
shp.df2$TOTAL <- as.numeric(shp.df2$TOTAL)
# 割合の計算
shp.df2$RATE <- shp.df2$TOTAL/shp.df2$JINKO
# エラー処理
shp.df2$RATE[is.infinite(shp.df2$RATE)] <- 0
shp.df2$RATE[is.na(shp.df2$RATE)] <- 0
#文字コード変換(SJIS > UTF8)
shp.df2$MOJI <- iconv(as.character(shp.df2$MOJI),from="SJIS",to="UTF8")
shp.text <- data.frame(shp.df2$X_CODE,shp.df2$Y_CODE,shp.df2$MOJI)
# 重複除去
shp.text <- unique(shp.text)
names(shp.text) <- c("lon","lat","name")
shp.text$name <- as.character(shp.text$name)
if(!file.exists("parking.xml")){
download.file(url="http://www3.city.sabae.fukui.jp/xml/parking/parking.xml",destfile="parking.xml")
}
#XMLの読み込み
doc<-XML::xmlParse("parking.xml",encoding="UTF8",useInternalNodes=TRUE)
#要素「parkinng1」を数える
len <- length(XML::getNodeSet(doc,"//parkinng1"))
#空のdata.frameを作る
dat <- data.frame()
#要素「parkinng1」の分繰り返す
for(i in 1:len){
# 要素値を取得する
label <- XML::xmlValue(XML::getNodeSet(doc,paste("//parkinng1[",i,"]/name"))[[1]])
latitude <- as.numeric(XML::xmlValue(XML::getNodeSet(doc,paste("//parkinng1[",i,"]/latitude"))[[1]]))
longitude <- as.numeric(XML::xmlValue(XML::getNodeSet(doc,paste("//parkinng1[",i,"]/longlatitude"))[[1]]))
size <- as.numeric(XML::xmlValue(XML::getNodeSet(doc,paste("//parkinng1[",i,"]/max"))[[1]]))
# 取得した要素値からdata.frameを作る
item <- data.frame(label,latitude,longitude,size)
dat <- rbind(dat,item)
}
# Googlemapを取得する
map <- ggmap::get_map(location=c(136.2302503,35.9633669),zoom=12,maptype="roadmap")
g <- ggmap::ggmap(map)
g <- g + ggplot2::geom_polygon(data = shp.df2, aes(x = long, y = lat, group = id,fill=JINKO),color="gray80",alpha=0.7)
g <- g + ggplot2::geom_point(data=dat,aes(x=longitude,y=latitude,size=size),color="blue")
g <- g+ ggplot2::scale_fill_gradientn(colours=rev(heat.colors(100)))
plot(g)
はじめに
オープンデータについてもう一度振り返ります． オープンデータとは「オープンなデータ」のことです． Open Definitionでは「誰もが自由に利用、再利用、再配布可能なデータ」と定義しています．
これまで演習では「オープンデータの評価指標「5 Star OPEN DATA（以下図参照）」に沿ったオープンデータの可視化に挑戦してきました．とくに三ツ星までの「オープンライセンス（Open License;OL）」「マシンリーダブル（Machine Readable;RE）」「オープンフォーマット（Open Format;OF）」の3要素が重要になります．
まず，そのデータがオープンであることを示すには，Creative Commons（[https://creativecommons.org/](https://creativecommons.org/） などのオープンライセンスを用いてオープンであることを明示する必要があります．
つぎに，データとして公開するのであれば， マシンリーダブル，すなわちコンピュータで処理しやすい形で提供する必要があります． PDFで公開されたデータより，Excelで公開されたデータの方が便利ですよね．
そして，データフォーマットもプロプライエタリ（データ形式等の仕様が非公開）なフォーマットよりは， オープンフォーマットの方が，処理可能性が広がり，誰もが利用できるデータとなるのです．
また，どんなデータでもオープンデータとなり得るのですが， その中でも価値あるデータと位置付けられているものがあります． G8オープンデータ憲章には「ハイバリューデータ・セット」として 利用価値の高いデータ種類がまとめられています．
| データカテゴリ | データセット例 |
|---|---|
| 企業 | 企業／事業者の登記情報 |
| 犯罪と司法 | 犯罪統計、安全 |
| 地球観測 | 気象／天候、農業、林業、水産業及び狩猟 |
| 教育 | 学校のリスト; 学校の業績、デジタルスキル |
| エネルギーと環境 | 汚染レベル、エネルギー消費 |
| 財政と契約 | 取引支出、賃貸契約、入札需要、将来入札、地方公共団体予算、国家予算（執行計画及び決算） |
| 地理空間 | 地勢、郵便番号、国の地図、地方（機関）の地図 |
| 国際開発 | 援助、食料安全保障、資源採掘、土地 |
| 政府の説明責任と民主主義 | 政府の窓口,選挙結果、法律及び規則、給与（給与水準）、福利厚生／贈与 |
| 健康 | 処方箋データ、実績データ |
| 科学と研究 | ゲノムデータ、研究及び教育活動、実験結果 |
| 統計 | 全国統計、国勢調査、社会基盤、財産、スキル |
| 社会的流動性と福祉 | 居住、健康保険及び失業手当給付金 |
| 交通とインフラ | 公共交通機関の時刻表、アクセスポイント、 ブロードバンドの普及度 |
ref. http://www.kantei.go.jp/jp/singi/it2/densi/dai4/sankou8.pdf
これまでこうした類のデータは専門家だけのものだったのですが， オープンデータやオープンサイエンスといったアイデアの浸透により，だれもが利用できるようになってきました．演習では「地球観測」「地理空間」「科学と研究」「統計」といったデータの可視化に挑戦してきましたが，なんとなくオープンデータの可能性については感じとってもらえたのではないでしょうか．
本日はRを使ったオープンデータの可視化の最終回となります． これまで以上に高度な内容が含まれていますが，挑戦する気概で演習に取り組んで下さい． 様々なデータの可視化については課題のヒントとしてご利用下さい．
Google Earthを使った可視化
これまでは ggplot2 や ggmap を使った可視化を試してきましたが， plotKML((http://plotkml.r-forge.r-project.org/)[http://plotkml.r-forge.r-project.org/])を用いると、Google Earthを用いた可視化が手軽にできます．
plotKML はその名のとおり，Google Earthで採用されているKML(Keyhole Markup Language)としてプロット結果を出力するライブラリです． 本演習では簡単な使い方を紹介します．興味がある方はチュートリアル(英語)，リファレンス（英語）等をご参照ください．
shapefileを使った可視化
さっそくですが，plotKMLの実行例です．統計GISからダウンロードしたshapefileを使って地域別人口を可視化します．
if(!require(maptools)){
install.packages("maptools")
library(maptools)
}
# KMLファイル用のライブラリ
if(!require(plotKML)){
install.packages("plotKML")
library(plotKML)
}
# shapefileを読み込む
shp<-readShapePoly("A002005212010DDSWC13212/h22ka13212.shp")
# 投影・座標系情報を設定する。
proj4string(shp) <- CRS("+proj=longlat +datum=WGS84")
# JINKOの値を高さに用いる
z<-shp@data$JINKO
z[z==0]<-1
# KML形式でplot
plotKML(shp,altitude=z)コードは単純ですが実行してみると， Google Earthが起動し，プロット結果を表示します．
加えて地名を表示してみます。
if(!require(maptools)){
install.packages("maptools")
library(maptools)
}
# KMLファイル用のライブラリ
if(!require(plotKML)){
install.packages("plotKML")
library(plotKML)
}
# shapefileを読み込む
shp<-readShapePoly("A002005212010DDSWC13212/h22ka13212.shp")
# 投影・座標系情報を設定する。
proj4string(shp) <- CRS("+proj=longlat +datum=WGS84")
# JINKOの値を高さに用いる
z<-shp@data$JINKO
z[z==0]<-1
#ラベル表示の準備
#文字コード変換(SJIS > UTF8)
shp@data$MOJI <- iconv(as.character(shp@data$MOJI),from="SJIS",to="UTF8")
# KML形式でplot
plotKML(shp["MOJI"],altitude=z)
散布図っぽく人口を可視化
今度はGoogle Earthで散布図もどきを描いてみましょう．
先ほどと同じように統計GISのshapefileを使います．
if(!require(maptools)){
install.packages("maptools")
library(maptools)
}
if(!require(plotKML)){
install.packages("plotKML")
library(plotKML)
}
shp<-readShapePoly("A002005212010DDSWC13212/h22ka13212.shp")
#必要な情報だけをもったdata.frameを作成する。
data <- data.frame(x=shp@data$X_CODE,y=shp@data$Y_CODE,id=shp@data$KEY_CODE,label=iconv(shp@data$MOJI,to="UTF8",from="SJIS"),jinko=shp@data$JINKO)
#座標を設定する。
coordinates(data) <- ~x+y
# 投影・座標系情報を設定する。
proj4string(data)<-CRS("+proj=longlat +datum=WGS84")
#KMLファイルを出力する。点（アイコン）のサイズは人口。
kml(obj=data,folder.name="",file.name="data.kml",shape="http://maps.google.com/mapfiles/kml/pal2/icon18.png",colour="#ff0000",alpha=1.0,size=jinko,labels=label)このように緯度経度情報を含んだdata.frameさえ用意すれば， ggplot2 と同じようにプロットすることができます．
shapeファイルに使えるアイコンセットは以下のサイトで確認することができます．
http://kml4earth.appspot.com/icons.html
data.frameを用意せず，kml関数に直接引数を渡すといった使い方も可能です．
if(!require(maptools)){
install.packages("maptools")
library(maptools)
}
if(!require(plotKML)){
install.packages("plotKML")
library(plotKML)
}
shp<-readShapePoly("A002005212010DDSWC13212/h22ka13212.shp")
proj4string(shp) <- CRS("+proj=longlat +datum=WGS84")
z<-shp@data$JINKO
z[z==0]<-1
kml(obj=shp,folder.name="",file.name="data2.kml",labels=iconv(shp@data$MOJI,to="UTF8",from="SJIS"),altitude=z,colour="#ff00ff",alpha=0.75,plot.lab=TRUE)
高齢者人口比率の可視化
最後に高齢者人口比率をplotKMLで可視化します．
if(!require(maptools)){
install.packages("maptools")
library(maptools)
}
if(!require(plotKML)){
install.packages("plotKML")
library(plotKML)
}
if(!require(plyr)){
install.packages("plyr")
library(plyr)
}
shp<-readShapePoly("A002005212010DDSWC13212/h22ka13212")
proj4string(shp) <- CRS("+proj=longlat +datum=WGS84")
csv<-read.csv("tblT000573C13212.txt",stringsAsFactors = FALSE,fileEncoding = "SJIS")
csv <- csv[-1,]
dat <- data.frame(KEY_CODE=csv$KEY_CODE,TOTAL=csv$T000573019)
shp@data <- join(shp@data,dat,by="KEY_CODE")
shp@data$TOTAL <- as.numeric(as.character(shp@data$TOTAL))
shp@data$RATE <- shp@data$TOTAL/shp@data$JINKO
# Inf(無限)となっていたものを0に置換(0除算結果の修正)
shp@data$RATE[is.infinite(shp@data$RATE)] <- 0
#NA(欠損値)を0に置換
shp@data$RATE[is.na(shp@data$RATE)] <- 0
#ラベル表示の準備
#文字コード変換(SJIS > UTF8)
shp@data$MOJI <- iconv(as.character(shp@data$MOJI),from="SJIS",to="UTF8")
#カラーパレットの準備
pal<-rev(heat.colors(100))
#確認
print(pal)
#カラーパレットを使って色を指定
kml(obj=shp,folder.name="",file.name="data3.kml",colour_scale=pal,colour=shp@data$RATE,labels=shp@data$MOJI,altitude=100,alpha=0.78,plot.lab=TRUE)ggplot2のときと可視化内容は同じですが、 Google Earthなのでインタラクティブに魅せることができるぶんインパクトがあります．
用途は限られますが，プレゼンテーション等の演出に使えるかもしれません．
HTML TABLEの可視化
ウェブ上にはTableタグを用いて公開されるデータも多く， Tableタグを処理できるようになると，より幅広くオープンデータが利用できるようになります．
RのXMLパッケージにはそのものずばりの readHTMLTable 関数があります． これを用いると，HTMLのTABLEタグをdata.frameとして読み込むことができます．
それでは具体例をもとにreadHTMLTable関数を使ってみましょう．
東京都の平均気温の可視化
気象庁のオープンデータを利用して東京都の平均気温の傾向を可視化します．
「過去の気象データの検索」にアクセスします．
http://www.data.jma.go.jp/obd/stats/etrn/index.php
「都道府県・地方を選択」をクリックします．
「東京」をクリックします．
「地点の選択」でも「東京」をクリックします．
「観測開始から月ごとの値を表示」をクリックします．
